% This function sets up an alternative state space representation 
% of the model solved up to fourth order where innovations only enter linearly.
% That is we transform the solution form
% z_t = f(x_t-1,u_t,sig)
% to 
% y_t   = g(x_t-1,u_t,sig)
% x_t   = h(x_t-1,u_t,sig)
% u_t   = eps_t
% where v_t = [x_t-1,u_t]
% By Martin M. Andreasen. 
function [g_ss,gv,gvv,gss,gvvv,gssv,gsss,gvvvv,gssvv,gsssv,gssss,...
    h_ss,hv,hvv,hss,hvvv,hssv,hsss,hvvvv,hssvv,hsssv,hssss,eta,ny,label_y,label_v,label_u] = ...
    StateSpace_form_LinearInov_4th(f_ss,fv,fvv,fss,fvvv,fssv,fsss,fvvvv,fssvv,fsssv,fssss,...
    sigma,nv,nu,nx,nz,dyn_vars,dyn_state_vars,dyn_shocks);

% We setup the function h. The dimensions are
h_ss    = zeros(nv,1); %the steady state level
hv      = zeros(nv,nv);
hvv     = zeros(nv,nv,nv);
hss     = zeros(nv,1);
hvvv    = zeros(nv,nv,nv,nv);
hssv    = zeros(nv,nv);
hsss    = zeros(nv,1);
hvvvv   = zeros(nv,nv,nv,nv,nv);
hssvv   = zeros(nv,nv,nv);
hsssv   = zeros(nv,nv);
hssss   = zeros(nv,1);
eta     = zeros(nv,nu);

% Setting eta
eta(nx+1:nv,1:nu) = sigma;

% This index keeps account of the variables from the f function which are put in the h function 
index_f = zeros(1,nz);

% We construct the h function
label_x = cell(nx,1);
for i=1:nx
    for j=1:nz
        if strcmp(cellstr(dyn_state_vars(i,:)),cellstr(dyn_vars(j,:))) == 1
            index_f(1,j)     = 1;
            label_x(i,:)     = dyn_state_vars(i,:);
            h_ss(i,1)        = f_ss(j,1);
            hv(i,:)          = fv(j,:);
            hvv(i,:,:)       = fvv(j,:,:);            
            hss(i,1)         = fss(j,1);
            hvvv(i,:,:,:)    = fvvv(j,:,:,:);            
            hssv(i,:)        = fssv(j,:);    
            hsss(i,1)        = fsss(j,1);    
            hvvvv(i,:,:,:,:) = fvvvv(j,:,:,:,:);            
            hssvv(i,:,:)     = fssvv(j,:,:);    
            hsssv(i,:)       = fsssv(j,:);
            hssss(i,1)       = fssss(j,1);
        end
    end
end
% The label for the shocks
label_u = dyn_shocks;
% The label for v
label_v = [label_x;label_u];

% We construct the g function. The dimensions are
ny    = nz-nx;
g_ss  = zeros(ny,1); %the steady state level
gv    = zeros(ny,nv);
gvv   = zeros(ny,nv,nv);
gss   = zeros(ny,1);
gvvv  = zeros(ny,nv,nv,nv);
gssv  = zeros(ny,nv);
gsss  = zeros(ny,1);   
gvvvv = zeros(ny,nv,nv,nv,nv);
gssvv = zeros(ny,nv,nv);
gsssv = zeros(ny,nv);   
gssss = zeros(ny,1);   
label_y = cell(ny,1);
i     = 0;
for j=1:nz
    if index_f(1,j) == 0
        i = i + 1;
        index_f(1,j)     = 1; %To indicate that we have used this position
        g_ss(i,1)        = f_ss(j,1);
        label_y(i,:)     = dyn_vars(j,:);
        gv(i,:)          = fv(j,:);
        gvv(i,:,:)       = fvv(j,:,:);
        gss(i,1)         = fss(j,1);
        gvvv(i,:,:,:)    = fvvv(j,:,:,:);
        gssv(i,:)        = fssv(j,:);
        gsss(i,1)        = fsss(j,1);
        gvvvv(i,:,:,:,:) = fvvvv(j,:,:,:,:);
        gssvv(i,:,:)     = fssvv(j,:,:);
        gsssv(i,:)       = fsssv(j,:);
        gssss(i,1)       = fssss(j,1);        
    end
end

end